#pragma once

#include "material.h"
#include "particle.h"
#include "point.h"
#include "type.h"
#include <numbers>
#include <tuple>

namespace PS
{

class ContactModel
{
public:
    std::tuple<Vector<2>, Vector<2>> normal_and_tangential_force(const Circle & p1,
                                                                 const Circle & p2,
                                                                 scalar delta_n,
                                                                 scalar delta_t,
                                                                 scalar e);

    scalar kn(const Circle & p1, const Circle & p2);

    scalar kt(const Circle & p1, const Circle & p2, double normal_overlap);

    scalar eta_n(const Circle & p1, const Circle & p2, scalar kn, scalar e);

    scalar eta_t(const Circle & p1, const Circle & p2, scalar kt, scalar e);

    scalar eta_r(const Circle & p1, const Circle & p2, scalar kr, scalar e);

    scalar collision_period_n(const Circle & p1, const Circle & p2, scalar kn, scalar e);

    scalar collision_period_t(const Circle & p1, const Circle & p2, scalar kt, scalar e);

    scalar collision_period_r(const Circle & p1, const Circle & p2, scalar kr, scalar e);

private:
    scalar equivalent_mass(scalar m1, scalar m2, scalar c);

    scalar eta(scalar equivalent_mass, scalar k, scalar e);

    scalar collision_period(scalar equivalent_mass, scalar k, scalar e);
};

} // namespace PS

namespace PS
{

inline scalar
ContactModel::kn(const Circle & p1, const Circle & p2)
{
    const auto miu1 = p1.material()->poisson_ratio();
    const auto miu2 = p2.material()->poisson_ratio();

    const auto e1 = p2.material()->elasticity_modulus();
    const auto e2 = p2.material()->elasticity_modulus();

    return 4.0 / 3.0 * 1.0 / ((1.0 - miu1 * miu1) / e1 + (1.0 - miu2 * miu2) / e2) /
           std::sqrt((p1.radius() + p2.radius()) / (p1.radius() * p2.radius()));
}

inline scalar
ContactModel::kt(const Circle & p1, const Circle & p2, double normal_overlap)
{
    const auto miu1 = p1.material()->poisson_ratio();
    const auto miu2 = p2.material()->poisson_ratio();

    const auto g1 = p2.material()->shear_modulus();
    const auto g2 = p2.material()->shear_modulus();

    return 8.0 * std::sqrt(normal_overlap) * ((1.0 - miu1 * miu1) / g1 + (1.0 - miu2 * miu2) / g2) /
           std::sqrt((p1.radius() + p2.radius()) / (p1.radius() * p2.radius()));
}

inline scalar
ContactModel::equivalent_mass(scalar m1, scalar m2, scalar c)
{
    return c * m1 * m2 / (m1 + m2);
}

inline scalar
ContactModel::eta(scalar equivalent_mass, scalar k, scalar e)
{
    const auto c = e == 0 ? 1.0
                          : std::log(1.0 / e) / std::sqrt(std::pow(std::numbers::pi, 2.0) +
                                                          std::pow(std::log(1.0 / e), 2.0));
    return 2.0 * std::sqrt(equivalent_mass * k) * c;
}

inline scalar
ContactModel::eta_n(const Circle & p1, const Circle & p2, scalar kn, scalar e)
{
    return eta(equivalent_mass(p1.mass(), p2.mass(), 1.0), kn, e);
}

inline scalar
ContactModel::eta_t(const Circle & p1, const Circle & p2, scalar kt, scalar e)
{
    return eta(equivalent_mass(p1.mass(), p2.mass(), 2.0 / 7.0), kt, e);
}

inline scalar
ContactModel::eta_r(const Circle & p1, const Circle & p2, scalar kr, scalar e)
{
    return eta(equivalent_mass(p1.inertia(), p2.inertia(), 1.0), kr, e);
}

inline scalar
ContactModel::collision_period(scalar equivalent_mass, scalar k, scalar e)
{
    const auto c =
        e == 0 ? 1.0 : std::sqrt(std::pow(std::numbers::pi, 2.0) + std::pow(std::log(e), 2.0));
    return 2.0 * std::sqrt(equivalent_mass / k) * c;
}

inline scalar
ContactModel::collision_period_n(const Circle & p1, const Circle & p2, scalar kn, scalar e)
{
    return collision_period(equivalent_mass(p1.mass(), p2.mass(), 1.0), kn, e);
}

inline scalar
ContactModel::collision_period_t(const Circle & p1, const Circle & p2, scalar kt, scalar e)
{
    return collision_period(equivalent_mass(p1.mass(), p2.mass(), 2.0 / 7.0), kt, e);
}

inline scalar
ContactModel::collision_period_r(const Circle & p1, const Circle & p2, scalar kr, scalar e)
{
    return collision_period(equivalent_mass(p1.inertia(), p2.inertia(), 1.0), kr, e);
}

inline std::tuple<Vector<2>, Vector<2>>
ContactModel::normal_and_tangential_force(const Circle & p1,
                                          const Circle & p2,
                                          scalar delta_n,
                                          scalar delta_t,
                                          scalar e)
{
    const auto k_n = kn(p1, p2);
    const auto k_t = kt(p1, p2, delta_n);
    
    const auto etan = eta_n(p1, p2, k_n, e);
    const auto etat = eta_t(p1, p2, k_t, e);

    const auto vref = p1.velocity() - p2.velocity();
    auto dirn = (p1.position() - p2.position()).normalize();
    auto dirt = dirn.rotate(90.0);
    auto fn = (-delta_n * k_n - etan * vref.inner_product(dirn)) * dirn;
    auto ft = (-delta_t * k_t - etat * vref.inner_product(dirt)) * dirt;

    return {fn, ft};
}

} // namespace PS
